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We examine the effectiveness of assuming an equal probability for states far from equilibrium. 
For this aim, we propose a method to construct a master equation for extensive variables describ- 
ing non-stationary nonequilibrium dynamics. The key point of the method is the assumption that 
transient states are equivalent to the equilibrium state that has the same extensive variables, i.e., 
an equal probability holds for microscopic states in nonequilibrium. We demonstrate an applica- 
tion of this method to the critical relaxation of the two-dimensional Potts model by Monte Carlo 
simulations. While the one-variable description, which is adequate for equilibrium, yields relaxation 
dynamics that are very fast, redundant two- variable description well reproduces the true dynamics 
quantitatively. These results suggest that some class of the nonequilibrium state can be described 
with a small extension of degrees of freedom, which may lead to an alternative way to understand 
nonequilibrium phenomena. 
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I. INTRODUCTION 

An surprising property of equilibrium states is that 
they can be identified by a small number of extensive 
variables (EVs) even though the system has an almost 
infinite number of degrees of freedom (DOFs). When EV 
values are given, the equilibrium system almost surely ex- 
hibits the most typical state, i.e., the equilibrium state is 
uniquely determined by macroscopic observation. This is 
the basis of the equilibrium statistical mechanics as the 
law of large numbers and an equal probability for micro- 
scopic states. On the other hand, nonequilibrium states 
generally need more information to be identified, and a 
history dependence (non-Markov process) appears if we 
label the state in a manner similar to that used to label 
the equilibrium states [1]. To treat the dynamics of a 
system as a Markov process, one has to deal with more 
DOFs, equal to the total number of the system in the 
worst case. However, it is natural to expect that there 
will be many situations in which a small extension of the 
DOFs will provide a sufficiently accurate description of 
nonequilibrium phenomena. In this paper, we propose 
a method to describe the non-stationary nonequilibrium 
dynamics by employing a few more EVs that are redun- 
dant for equilibrium states. We start from a macroscopic 
description with a small number of DOFs and approach 
the microscopic description by increasing degrees. This 
approach proceeds in the manner opposite to the usual 
approach, which start from microscopic description and 
reduces the DOFs although approaching the same point. 
Thus, this approach may be a new method to connect 
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the microscopic and macroscopic worlds and add a new 
perspective to nonequilibrium statistical physics. 

If out-of-cquilibrium dynamics can be described by the 
information of equilibrium states, like the successful lin- 
ear response theory [2], it can be very useful for our un- 
derstanding. In the vicinity of the equilibrium, we have 
the phenomenological theory stating that the thermody- 
namic force that restores the system to an equilibrium 
is proportional to the free energy gradient, which is a 
function of EVs [3-5]. A quadratic well of the free en- 
ergy is often assumed, which leads to a linear response to 
the deviation from its lowest point, i.e., the equilibrium 
point. This theory is useful for a weak nonequilibrium, 
but cannot be applied for a strong nonequilibrium. 

The projection operator method [1, 6, 7] is typically 
used to extract the dynamics of slow-changing variables, 
i.e., integrating out the rapidly changing variables. In 
the formalism by Mori [1], the projection leads to a gen- 
eralized Langevin equation as an equation of motion for 
the remaining variables. The reduction of the number of 
DOFs can be seen as a local equilibrium approach. We 
take a similar approach; specifically, we develop a method 
to describe stochastic dynamics. 

Consider the case of the relaxation dynamics far from 
equilibrium, where intensive variables such as tempera- 
ture are controlled. It is difficult to relate a general tran- 
sient state to the equilibrium state that corresponds to 
the external condition. However, a transient state labeled 
by a set of EVs, such as internal energy, can be related to 
another equilibrium state having the same energy. It is a 
microcanonical ensemble that is equivalent to the canon- 
ical ensemble at a temperature different from that of 
the environment. As mentioned before, nonequilibrium 
states need a higher number of EVs than equilibrium 
states. Here we call such states "extended equilibrium 
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states" . The equivalence of the transient state to the ex- 
tended equilibrium state means that an equal probability 
holds for the noncquilibrium state. It is worth consider- 
ing the possibility of describing nonequilibrium dynamics 
by using only information from the extended equilibrium 
state. If such description is possible, it will be a useful 
clue for understanding noncquilibrium phenomena based 
on the well-understood equilibrium property. 

This paper is organized as follows. In the next section, 
we propose a method to construct a master equation by 
using only the equilibrium and macroscopic information. 
In Sees. Ill, we analyze the all-to-all coupling Ising model 
as the simplest example. In Sec. IV and V, we apply our 
approach to the three-state Potts model on the square 
lattice with one- and two-variable descriptions, respec- 
tively. Section VI presents the discussion and perspec- 
tive. 



II. METHOD 
A. Master equation 

We consider a dynamics represented by a master equa- 
tion in terms of a microscopic state variable X = 
(CTi,CT 2 ,---,CTA r ) as 



side becomes 

£ w(b;X'\X)p(X,t)6 A A 



X,X' 



= E E %,iMb;X'\x)s IxI p(x,t) 

A' X,X' 

= Y^W{b;A'\A)P{A,t). (4) 
A' 

To obtain the last line, we use Eq. (3) and put 

W (b;l'\l) ee -L- £ S A A ,w(b;X'\X)6 A A 
9(A) X*, 



= (Y,hf,AMb-,X'\X)) , 



(5) 



X' 



where the bracket indicates the unweighted average with 
A fixed as 



(6) 



By performing the same operation to the first term in the 
right hand side of Eq. (1), we obtain the following form 
closed for P(A, t) as 



dp(X,t) 
dt 



£ [w(b;X\X')p(X',t) - w(b;X'\X)p(X,q}) ^AlH = £ [ W (b;A\A')P(A',t) - W(b; A'\A)P(A, t\]) 



where p{X,i) is the probability distribution function 
(PDF) of X at time t, and w(b;X'\X) is the transition 
probability per unit time from X to X' with given b. 
Here b — (b±, 62, • • ■ b n ) is a set of intensive variables, 
such as inverse temperature, magnetic field, or chemical 
potential, introduced as an external control parameter. 
Our aim is to derive the corresponding master equation 
in terms of the EVs A = (Ax, A2, ■ ■ ■ , A n ) that are con- 
jugate to 6, such as internal energy, magnetization, or 
number of particles. A is a function of X, which we de- 
note as A = Ag; thus, the PDF of A is given by 



p(l,t) = J2p(x,t)s A A . 



(2) 



The main assumption of the present formulation is the 
equal probability for microstates for any given A and t 



p(X,t) = P(Ag,t)/g(A), 



(3) 



where g(A) is the number of states having A as g{A) = 
Y,x 5 A x A- B y multiplying Eq. (1) by S A A and tak- 
ing the summation over X, the left hand side yields 
dP(A,t)/dt. This represents a projection from the X- 
space to the A-space. The second term on the right hand 



It is often the case that the transition probability w is 
expressed as 

w(b;X'\X) =w(b;Ajf, - A^)D(X'\X), (8) 

where the function D(X'\X) equals unity if a direct (one- 
step) transition path exists between X and X', and is 
zero otherwise. For example, w(b; AA) is proportional to 

min[l,e~ &A ^] with AA = Ag, — A% in the Metropolis 
algorithm. The expression Eq. (8) reflects the fact that 
the change of A is decomposed to the local change by the 
short-range interaction, which is not explicitly dependent 
on the accumulated value A itself. In such a case, W is 
expressed as 

W(b;A + AA\A) = fi(A; AA)w(b; AA) (9) 
with ^{A-M) = {^5 AtiA+AA D{X'\X)\ .(10) 

\ X> I A 

Here fi(A; AA) is regarded as the (un-normalizcd) PDF 
of the possible path, causing AA for a given A. Note that 
(i itself is not a dynamical variable but a state quantity 
that is statistically evaluated in the extended equilibrium 
ensemble mentioned before. Thus, the information of 
the spatial configuration is reduced to a small number of 
quantities: the components of AA. 
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Since the general nonequilibrium state requires a 
macroscopic number of DOFs, the present master equa- 
tion to resemble a Markov process with finite DOFs is an 
approximated one. However, it is possible to improve its 
accuracy in a systematic manner by increasing the num- 
ber of EVs. In our strategy, [i is estimated from the equi- 
librium ensemble: an equal probability with fixed EVs A. 
Although the factorization of Eq. (9) is practically very 
useful and naturally derived for standard Monte-Carlo 
(MC) dynamics, it is not clear whether it can be used 
for any system. If it cannot be used, we have to estimate 
W(b; A + AA\A) for every pair of b and A as in Eq. (5). 

One may think that the dynamics with the transition 
probability of Eq. (9) with w(b; AA) oc min[l, e _& A4 ] is 
similar to that with 

W(b; A'\A) oc min[l, e -[F$*)-F&i$]] (11) 

where F{b; A) = b ■ A — In g(A) is the dimensionlcss free 
energy, and g{A) is the number of states having A. Such 
an approach was examined in Refs. [8-10]. Although 
both of Eqs. (9) and Eq. (11) lead to the equilibrium 

PDF P cq oc e~ Fi < b ' A ^ as a stationary distribution, the 
time evolutions are different. Our strategy can construct 
the equation of motion for the EVs from first-principles, 
as shown below. The resultant equation depends on the 
function form of w(b, AA) and does not includes a fitting 
parameter, such as a dissipation coefficient. 

B. Equation of motion for the EVs 



Hereafter, we only treat such most probable dynamics. 

If the rotation (d/dA) x V equals zero, we can define 
a potential function F ncq (b; A), which can be regarded as 
nonequilibrium free energy, as 

V(b;A) = -^F ncq (b;A), (14) 
oA 

which depends on the details of the dynamics represented 
by w(b; AA). This is in contrast to the dynamics corre- 
sponding to Eq. (11), which is governed by the equilib- 
rium free energy. 

C. Equilibrium ensemble with fixed EVs 

To construct the master equation, we need to calcu- 
late \i for various A values as a statistical distribution 
in the extended equilibrium ensemble. We calculate it 
by MC simulations. It is not efficient to perform an in- 
dependent simulation for each A. Instead, we perform 
the Wang-Landau sampling [11, 12] to obtain a flat his- 
togram of A. The Wang-Landau sampling enables us 
to cover a wide range of A uniformly by a single run. 
In addition, a joint density (number) of states g(A) is 
obtained as a byproduct, which leads to the microcanon- 
ical entropy hig(A). The generalized dimensionlcss free 
energy is given by F(b; A) = b ■ A — \ng(A), and its ex- 
tremal condition, dF/dA — 0, gives an equilibrium rela- 
tion, b = d\ng(A)/dl [13]. 



Let n be the dimension of b and A. At equilibrium, 
A is a unique function of b except for the multistable 
points. [On the other hand, b is a unique function of 
A even in coexisting phases.] The equilibrium state of 
a given b is represented by a point in the n-dimensional 
space, which we call the A-space. On the other hand, 
the trajectory of nonequilibrium dynamics is given by a 
curved line. Generally, each point of the trajectory, i.e., 
a transient state, is not related to the equilibrium state 
with b of the environment. 

We can define the velocity field V = (Va 1 , • • • , Va„ ) as 



V(b; A) = AAW(b; A + AA) 



(12) 



AA 



which is used to express the probability current as 
3{b;A,t) = P(A,t)V(b;A). Note that V is independent 
of time. The dynamics is deterministic in the thermo- 
dynamic limit, where the fluctuation of the EVs can be 
ignored in comparison with the mean values. The time 
evolution of (A) (t) = J dAAP(A, t) can then be given as 
a solution of the equation of motion, 



d(A) 
dt 



V(b; (A)). 



(13) 



III. THE ALL-TO-ALL COUPLING ISING 
MODEL 

Here we demonstrate how the above method is applied 
to an analytically solvable model. Consider a two-state 
Potts model, i.e., the Ising model containing N spins with 
all-to-all coupling. The energy function is 



E 



1 



V 51 



±[N* + (N ~ N Q ) 2 }. (15) 



The spin variables {<Ji} take the value or 1. The number 
of spins at state is defined as 



N = J2 S ^o 



Nn . 



(16) 



Since the present system does not have a meaningful spa- 
tial structure, the state can be accurately identified only 
by Nq. From the ./^-symmetry, it is sufficient to consider 
No > N/2. Therefore, the magnetization and energy 
cannot change independently but have a one-to-one re- 
lation, unlike infinite-dimensional systems. By utilizing 
this property, we use the notation as if Nq were conjugate 
to inverse temperature /3 in the following; E and AE as 
the arguments of w and fi are replaced with N and A No, 
respectively 
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The number of states g(No) is a function of Nq as 

Nl 



g(N ) 



(17) 



N \(N-N )V 
and the dimensionlcss free energy per site is written as 



^ « -fi[nl + (1 - n ) 2 } 



hv 



1 



1 -n 



no In 



no 



1 - n 



(18) 



= In 2 - | + i(l - 0)m 2 + ^m 4 + 0( ? n 6 ), (19) 

where = is an inverse temperature and m = 

2tiq — 1 is the order parameter. It is well known that 
this is equivalent to the $ 4 free energy, i.e., the Landau 
model, which has a single minimum at m = for < 
C = 1 and has two symmetric minima at ±m(/3) with 
< m(0) < 1 for > P c ; m{0) = ±^Z{0- 1) for 
m <C 1. 

We consider the Metropolis dynamics; randomly pick- 
ing up one spin and changing the state to 1 if it is and 
vice versa. Therefore, the updated PDF is 



u{N ; +1) = 1 — no and u(yV ; -1) = n . 
The spin flip is accepted with probability rate 

N 



w(b; AA) 



min[l, e 



-b-AA-i 



(20) 



(21) 



where q — 2 for the Ising model. For No > N/2, 

w(P; +1) - N, w(0; -1) = Ne~ pAE (22) 
with AE = 2{2n Q -l) + 0(N- 1 ). (23) 

The velocity can now be written as 

vn„ = V No /N 

= fx(N ;+l)w(p; +1) - u(N ; -l)w(0; -1) 

= l-n [l + e^ 2 ' 3(2,lo - 1) ]. (24) 

The time evolution is formally solved as 

dn 



t = N 



no(0) VN (P,n )' 



(25) 



At the critical point = 1, we have Ujy = — m 3 + 0(m 5 ), 
which leads to m cx t~ x / 2 for m «l. The time evolution 
of energy per site behaves as (E c — E)/N rts m 2 /2 oc 
where E c — N/2 is the equilibrium energy at /3 C . 

We now consider the relation between the above dy- 
namics and the free-energy landscape. The free-energy 
gradient is written as 



dF 

dm 



2/3(2n - 1) 



In- 



fo 



1 - no 



(26) 



We have -dF/dN = -4m 3 /3 + 0(m 5 ) at = 1. Thus, 
the phcnomcnological equation of motion, dNo/dt = 



master eq. 




FIG. 1. Time evolutions of order parameter m = 2no — 1 
(multiplied by t 1 ^ 2 ) in the all-to-all coupling Ising model. The 
decimal logarithm is shown on the horizontal axis. The result 
of the exact solution is denoted by "master eq." and that 
obtained form the free-energy gradient is denoted by "free 
energy" . While the long-time behaviors of the two results are 
same as m oc t -1 ' 2 , the short-time dynamics are significantly 
different. 



—dF/dN , also yields m oc i~ 1//2 . The short time dy- 
namics is, however, distinctly different from the exact so- 
lution, Eq. (25). Especially, dF/dN diverges for n = 
and 1, as against the exact solution. Figure 1 shows 
the results of the numerical integrations of Eqs. (25) and 
(26) with an initial condition no = 0.95 at t = 0. It is 
observed that the approximated dynamics with the equi- 
librium free energy is valid only in the vicinity of the 
destination state, i.e., the equilibrium. 



IV. ONE- VARIABLE DESCRIPTION 

Next, we investigate a finite-dimensional system: the 
ferromagnetic 3-state Potts model on the square lattice, 
which exhibits a second-order transition at the critical 
point C = 111(1 + ^3) « 1.00505. Here we did not choose 
the Ising model (q = 2) to avoid an unusual critical be- 
havior, the logarithmic divergence of specific heat. 

We start with the minimal case, a description with one 
EV. For a system with fixed temperature, the most rele- 
vant EV is internal energy E, which is an EV conjugate 
to temperature. If we do not use it, the transition rate, 
which is a function of 0E, cannot be determined. The 
energy E is defined as 



E = 2N- <W, 

(i,j)Gn.n. 



(27) 



where the summation is taken over all the nearest- 
neighbor pairs. Each spin Gi takes a value: 0, 1, or 2. 
The energy change of a single spin flip is in the range of 
—4 < AE < 4. We performed the Wang-Landau sam- 
pling simulations [12] for samples with size N = LxL and 
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FIG. 2. (color online) (a) Energy dependence of velocity for 
the one- variable description. Colors denote temperatures, (b) 
Velocity observed in the KMC. Colors denote temperatures. 
Non-monotonic behaviors are observed and finite-size depen- 
dence appears near the zero-crossing points of Ve- 



L = 64 — 1024 by imposing a periodic boundary condi- 
tion. We averaged the PDF (J,(E; AE) over four indepen- 
dent runs. In the estimation of fj,, we eliminated the up- 
date that the spin state remains the same. We also per- 
formed corresponding kinetic Monte Carlo (KMC) sim- 
ulations [14, 15] as a non-approximated dynamics and 
compared the results with those obtained by the master 
equation. One MC step involves the following process N 
times: picking up a spin randomly and flipping it to one 
of the q—1 states with probability min[l, e~ b ' ], Thus, 
the transition rate of AA ^ is written as Eq. (21). 

Figure 2(a) shows the velocities Ve(E) at several tem- 
peratures, which were obtained by the master equation 
with the one-variable description. Each zero-crossing 
point of Ve is related to the equilibrium state at each 
temperature. From these velocities, the time evolutions 
of the energies are obtained by the numerical integration 
like Eq. (25) and are plotted in Fig. 3(a) along with the 
results of the KMC simulations. While the results ob- 
tained by the two methods have similar the equilibrium 
energy, the relaxations obtained by the master equation 
are significantly faster than those obtained by the KMC 



FIG. 3. (color online) (a) Time evolution of energy. Colors 
denote temperatures. The decimal logarithm is shown on the 
horizontal axis. The results of the one- variable description are 
denoted by lines and those of KMC simulations are denoted 
by symbols. While the equilibrium values are the same, the 
short-time dynamics are significantly different between the 
two methods, (b) Relaxation of energy at the critical tem- 
perature. The decimal logarithm is shown on the both axes. 
Colors denote temperatures. Three groups are shown: the 
KMC (symbol), the one- variable description (dashed line), 
and the two- variable description (solid line). 



method. 

Figure 2(b) shows the velocities calculated by the PDF 
/i obtained by nonequilibrium sampling through KMC 
simulations. These velocities reproduce the time evolu- 
tion of the energy obtained by the KMC method. Non- 
monotonic behaviors observed just before Ve crosses zero 
level arises as a result of the crossover from nonequilib- 
rium sampling to equilibrium sampling in the long time 
regime, and the latter sampling coincides with the equi- 
librium velocity in Fig. 2(a). When the system size is 
increased, this equilibrium region becomes smaller be- 
cause the equilibrium fluctuation of E/N disappears for 
N — > oo. In the thermodynamic limit, the velocities 
of equilibrium and nonequilibrium sampling are different 
for all E except for the equilibrium point. This means 
that the typical spin configuration of a given E, which 
is represented by /i, is different between the two states. 
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x = |E-E |/N 

c 

FIG. 4. (color online) Confirmation of the power-law at 
P = /3 C . The microcanonical velocity (a) and microcanonical 
temperature (b). It is shown that \Ve\ and |/3 — /3 C | converge 
to finite values by multiplying x 1 " 3 , where x = \E — E c \/N, 
and taking the limit x — > 0. This means the two quantities 
are proportional to \E — E c \ ■ 



Thus, a one- variable is insufficient to describe this type of 
dynamics although it yields the correct equilibrium state. 

As mentioned above, the dynamics obtained by the 
one-variable description fail even in the vicinity of the 
equilibrium point. Here we look at it more closely. Fig- 
ure 4(a) shows the relation between Ve and E at f3 c , 
which indicates a power-law, \Ve\ oc \E c — E] 1 ' 5 , in the 
limit: E — > E c . This leads to E c - E oc t~ 2 - although 
it is not confirmed in the short-time relaxation as shown 
in Fig. 3(b) [dashed lines]. The exponent 2.0 is distinctly 
different from that obtained by the KMC method, 0.367 
[16, 17]. 

Next, we examine the relation between the velocity 
and the free-energy gradient, —dF((3 c ; E)/dE = /3(E) — 
p c . The microcanonical temperature (3(E) is defined as 
d\ng(E)/dE. Figure 4(b) shows a power law singular- 
ity 13(E) - (3 C oc (E - -Ee) 1 ^ 1 ""), where a is the criti- 
cal exponent of the specific heat. The exponent is con- 
sistent with the known value 3/2 (a = 1/3). There- 
fore, the phcnomcnological dynamics using free energy 
also fail to reproduce the true dynamics in a similar 
manner to that using Ve- We can write the relation 
Ve(E) oc —dF(j3 c ;E)/dE since both the velocity and 
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FIG. 5. (color online) (a) Time evolution of the internal en- 
ergy. The results using the two-variable description (lines) 
and those obtained by the KMC method (symbols) are in 
good agreement, (b) Time evolution of the order parameter. 
The line with the slope P/zv = 0.0602 [17] is also shown. The 
two- variable description reproduces the critical behavior well. 



the gradient are proportional to (E — Ec) 1 - 5 . Therefore, 
the noncquilibrium free energy defined in Eq. (14) is es- 
sentially equivalent to the equilibrium energy in the one- 
variable description. This will be discussed more closely 
after seeing the two-variable description in the next sec- 
tion. 



V. TWO- VARIABLE DESCRIPTION 

Next, we employ one more EV, iVo, that is the number 
of spins taking the 0-state, which is related to the order 
parameter m = (Nq/N — l/q)/(l — 1/q). The extensive 
and intensive variables are given by A — (E,Nq) and 
b = ((3, —h), respectively. Here h is the magnetic field 
coupled with the 0-state as —hNo- In the following, we 
set h = 0. The update of No with a single spin flip is 
restricted to — 1 < ANq < 1. For the two dimensional A- 
space, we calculate n(E, M; AE, AN ) whose total num- 
ber of components is proportional to iV 2 x 9 x 3. This 
requires such a large memory space during the Wang- 
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E/N 

FIG. 6. (color online) The trajectory in (E x iVo)-space 
of /3 e xt = 1.010 > /3 C with various initial conditions. The 
fixed point is indicated by the open circle at (E,No)/N = 
(0.370, 0.805). The dashed line denotes h(E, No) = 0. 

Landau sampling that we divided E and N into bins of 
width Due to this coarscgraining, the integra- 

tion with the initial condition E = and N = N was 
difficult because the trajectory runs near the upper edge 
of the A-space, iVo = 1 — E/4, above which there is no 
state. As an initial condition, we used the values of E 
and No obtained by the KMC simulations at t = 4 for 
each temperature. 

We integrated the equation of motion, Eq. (13), with 
the 2nd order Rungc-Kutta method with a discrete time 
step of 0.01. Because we had V values only on the lattice 
points, we interpolated the value at the off- lattice points 
as V(E,N ) = Co + C\E + c 2 7V , where the coefficients 
were determined by the three nearest lattice points. 

The time evolution of E is shown in Fig. 5(a). The 
results of two different system sizes L = 64 and 128 are 
shown as dashed and solid lines, respectively. In a regime 
where the two lines are identical, the behaviors can be 
considered to be in the thermodynamic limit. The re- 
sults of the two-variable description show a reasonably 
better agreement with those of the KMC method than 
the one- variable description result do (the relaxations are 
still slightly faster than those of the KMC). The critical 
behavior is shown in Fig 3(b) as solid lines. We found 
that the relaxations obtained by the master equation 
agree reasonably well with those of the KMC method. 
The time evolutions of the order parameter are shown in 
Fig. 5(b). The relaxations obtained by the master equa- 
tion and KMC are in good agreement as is the internal 
energy. 



A. Trajectory 

Figure 6 shows the flow diagram in the (E x 7Vo)-spacc 
at j3 = 1.010 > j3 c . The flow line is drawn by trajecto- 
ries that are obtained by integrating V for various initial 
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FIG. 7. (color online) (a) Trajectories for the system L — 128. 
The results obtained by the integration of the master equation 
are denoted by lines and those obtained by the KMC method 
are denoted by symbols. The results of the two-variable de- 
scription reproduce the true dynamics significantly well. The 
contour line h(E,No) = (dashed line) is also shown, (b) 
The contour line of /3(A). 



conditions. The trajectories rapidly join the unique one 
flowing into the fixed point. Because of the rotation-free 
velocity field, the nonequilibrium free energy Eq. (14) 
can be defined, and it differs from the equilibrium one. 
In Fig. 7(a), trajectories obtained by the master equation 
and KMC are plotted and show good agreement for each 
temperature. 

The line h(E,N ) = -din g(E, N )/dN = is shown 
in Fig. 6. This line represents the state equation without 
the magnetic field. [The left below region of the line is 
the coexisting phase of domains ordered to states 0, 1, 
and 2. To be precise, h(E,No) also equals zero in this 
region in the thermodynamic limit.] Each point on this 
line is related to the most probable state in the fixed E 
ensemble, i.e., the microcanonical ensemble, without the 
magnetic field. Therefore, in the one-variable descrip- 
tion, the transient state is confined to this equilibrium 
line independently of (3. On the other hand, the trajec- 
tory of the disordering process in the two-variable de- 
scription runs above and to the right of the equilibrium 
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line and stops when touching the line. By adding the 
second variable iVo, the system can take a path away 
from the equilibrium line depending on p. This is the 
reason for the drastic improvement in accuracy from the 
one- variable description. 

B. Comparison with the free energy valley 

While the contour line h(E,No) = is related to 
the equilibrium line in the fixed-i? ensemble for h = 0, 
0(E,Nq) — Pcxt is related to the equilibrium line for the 
fixed- Nq ensemble for a given /3 e xt- The latter gives a 
trajectory of the quasi-equilibrium dynamics when the 
relaxation of the magnetization is much slower than that 
of the energy. Figure 7(b) shows the contour lines of 
p. The contour lines are similar to the trajectory of the 
isothermal relaxation obtained by the KMC method al- 
though the agreement is slightly less than that of the mas- 
ter equation results [Fig. 7(a)] far from equilibrium. The 
quasi-equilibrium dynamics cannot reproduce the rapid 
change in Nq, but that has an almost negligible effect in 
the present case. This is a useful approximation for the 
trajectory; note, however, that the contour lines of the 
temperature provide no information on the time evolu- 
tion. 



VI. SUMMARY AND DISCUSSION 

In this paper, we propose a method to analyze the 
non-stationary nonequilibrium state by using a master 
equation for extensive (macroscopic) variables and show 
some applications of the method. The method uses the 
quantity n(A; AA) calculated by equilibrium ensemble for 
fixed EVs, which means an equal probability is assumed 
for microscopic states in nonequilibrium. 

For the all-to-all coupling model, the state can be iden- 
tified accurately only by the single variable Nq since the 
system lacks a spatial structure. Thus, there is no differ- 
ence between equilibrium and nonequilibrium states in 
this case, and the one- variable description gives rigorous 
dynamics. For the two-dimensional 3-state Potts model, 
while the one- variable description is far from adequate to 
resemble the true dynamics, the two- variable description 
is sufficient. This improvement can be attributed to the 
fact that we have a unique trajectory for any external 
conditions (3 in the one-dimensional ^4-space correspond- 
ing to E; this trajectory is equivalent to the equilibrium 
line in two dimensions (this is the reason why strange 
inflection points, which corresponds to the singularity at 
P = p c , are observed for off-critical temperatures as in 
Fig. 2(a).) It has been suggested that transient states in 
finite-dimensional systems have distinctly different kinds 
of spatial fluctuations depending on temperature even for 
the same internal energy. Specifically, the trajectory runs 
in the region where h(E, N ) > 0. For the same E, N on 
the trajectory is larger for higher temperatures. This is 



because the minority spins with number (N — No) have a 
lesser tendency to aggregate at higher temperatures and 
thus acquire E with smaller (N — No) from the surface 
energy. 

To investigate equilibrium critical properties from 
nonequilibrium relaxation at the critical point, the so- 
called nonequilibrium relaxation analysis [18] can be 
used. In the relaxation, a nonequilibrium correlation 
length grows with time in a power law, which is sim- 
ilar to the growth of the equilibrium correlation length 
when approaching the critical point by tuning an external 
parameter. The critical exponents are estimated using 
this similarity [19]. The two paths on which the growth 
of the nonequilibrium and equilibrium correlation length 
are observed correspond to the trajectory of P = p c and 
the equilibrium line h(E,No) = 0, respectively. Both 
paths approach the critical point (E,Nq) = 7V(e c ,l/3), 
where e c = 1 — 1/V3, and the paths are expressed as 
{N /N - l/Zf/P = C{e c - £/A0 1/(1 ~ Q) - However, the 
proportionality constant C is different between the two, 
which makes the two paths macroscopically distinguish- 
able. Interestingly, while the relaxation dynamics along 
them are different, the critical exponents, such as a and 
P, are the same. 

For the two-dimensional 3-state Potts model, the sec- 
ond variable, the order parameter, has special impor- 
tance because its relaxation is much slower than that 
of the energy. Thus, it has more impact on the total dy- 
namics. Therefore, the quasi-equilibrium dynamics with 
respect to the order parameter is a good approximation. 
Indeed, the contour line of the temperature, i.e., the ar- 
ray of equilibrium points for various temperatures, gives 
a good analog of the relaxation trajectory it does not 
provide information on the time evolution. One would 
think that the one-variable description with No would 
give good results comparable to those of the two- variable 
description with E and Nq. However, in our formulation, 
the EV conjugate to the non-zero intensive variable must 
be employed to determine w(b; AA). 

It is possible to improve the accuracy of our method 
by employing more DOFs; however, the difficulty is to 
find appropriate variables independent of each other. 
Important criteria are the slowness of the relaxation 
and the capture of the spatial structure. One good 
candidate would be the structure factor, S*o(q) = 
I J2i <5o-,oe 4q ' ri | 2 /A^. The quantity N a is equivalent to 
\/NSo(0). By adding So from small q to larger ones, 
a finer structure of the space would be included stepwise. 
However, in practice, the addition of DOFs is very dif- 
ficult in the implementation adopted here because the 
dimensionality of the A-spacc increases. Depending on 
the problem, we might not need information of the wide 
range of the A-space. In addition, the calculation of \x is 
suitable for trivial parallel computation. 

One relaxation dynamics of a given intensive variable 
can be generated with a reasonably lower computational 
cost by KMC method than by the present method. Our 
method takes significant computational time to estimate 
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n(A; AA). Once this is obtained, however, one can easily 
generate time evolutions for various b. This is similar 
to the efficient calculation method of equilibrium expec- 
tation values by using the reweighting method with ex- 
tended ensembles [12]. This aspect is very useful when 
calculating the fluctuation of quantities, which requires 
many samples. However, our main aim is not to develop 
an efficient way to yield time evolution but to find a 
way to understand the system far from equilibrium in 
the framework of equilibrium statistics. The basic con- 
cept proposed here is not only for numerical analysis but 
for analytical treatment. It is possible to develop an ap- 
proach to noncquilibrium dynamics by utilizing the at- 
tained knowledge of equilibrium physics. The success in 
the present demonstration encourages us to continue in 
this direction. 

We have found that the present method works well for 
critical relaxation dynamics. The next interesting appli- 
cation is the dynamics taking first-order transitions [20] , 
such as liquid- vapor transition and protein folding [21]. 
It is often explained that the nucleation dynamics is gov- 



erned by the free-energy barrier determined by the com- 
petition between bulk and surface free energies. However, 
for a finite-dimensional system with short-range interac- 
tion, this picture is not applicable to the macroscopic free 
energy, which is a convex function of EVs [13] . Therefore, 
the addition of DOFs is essential for nucleation dynam- 
ics. 

The basic idea of the present method can be applied 
to a non-stationary state of an open system with cer- 
tain driving forces; a non-stationary state is described 
by a noncquilibrium stationary state. By replacing the 
word "equilibrium" with "stationary" , we obtain a par- 
allel framework for transient dynamics toward stationary 
states. 
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